%%  multipole scattering powerd=30um
clc; clear all; clf; close all;
c = 3e8; % %光速[m/s]
% n = 51; % % n为取的频率点数
p = 8814; % % CST中每一个频率提取的坐标点数
namePool = [100:2:216 217:227 228:2:300];
f = [1.1:0.002:1.216 1.217:0.001:1.227 1.228:0.002:1.3] * 1e12; % % Hz 整段频率
w = 2 * pi * f; % %角频率

%% 导入数据，保存成cell
for n = 1:length([100:2:216 217:227 228:2:300])

    str = ['.\SC\SC__', num2str(namePool(n)), '.txt'];
    importedDataName{n}=importdata(str);
    importedDataValue{n} = importedDataName{n}.data;
    kx{n} = complex(importedDataValue{n}(:, 4), importedDataValue{n}(:, 5));
    ky{n} = complex(importedDataValue{n}(:, 6), importedDataValue{n}(:, 7));
    kz{n} = complex(importedDataValue{n}(:, 8), importedDataValue{n}(:, 9));
    x = importedDataValue{1}(:, 1) * 1e-6; %x坐标
    y = importedDataValue{1}(:, 2) * 1e-6; %y坐标
    z = importedDataValue{1}(:, 3) * 1e-6; %z坐标
    s = importedDataValue{1}(:, 10) * 1e-12;
    Kxx(:, n) = [kx{n}(:, 1)];
    Kyy(:, n) = [ky{n}(:, 1)];
    Kzz(:, n) = [kz{n}(:, 1)];
end

xx = repmat(x, 1, n); % %x坐标重复扩充n列,[m]
yy = repmat(y, 1, n); % %y坐标重复扩充n列,[m]
zz = repmat(z, 1, n); % %z坐标重复扩充n列,[m]
S = repmat(s, 1, n); % %面积s重复扩充n列,[m^2]
r = [x, y, z]; % %坐标写成位矢,[m]
R = repmat(r, 1, n); % %位矢r重复扩充n组,[m]

%% 电偶极矩
Px = 1 ./ (i .* w) .* dot(S, Kxx);
Py = 1 ./ (i .* w) .* dot(S, Kyy);
Pz = 1 ./ (i .* w) .* dot(S, Kzz);

coefficient_P = ((2 .* w.^4) / (3 .* c^3));
Ipx = coefficient_P .* (abs(Px)).^2; % %电偶极矩辐射能量
Ipy = coefficient_P .* (abs(Py)).^2;
Ipz = coefficient_P .* (abs(Pz)).^2;
Ip = sqrt(abs(Ipx).^2 + abs(Ipy).^2 + abs(Ipz).^2);
clearvars Px Py Pz cofficient_P Ipx Ipy Ipz

%% 磁偶极矩
mx = (1 ./ (2 * c)) .* dot((conj(Kzz) .* yy - conj(Kyy) .* zz), S);
my = (1 ./ (2 * c)) .* dot((conj(Kxx) .* zz - conj(Kzz) .* xx), S);
mz = (1 ./ (2 * c)) .* dot((conj(Kyy) .* xx - conj(Kxx) .* yy), S);
% mx = (1 ./ (2 * c)) .* dot(((Kzz) .* yy - (Kyy) .* zz), S);
% my = (1 ./ (2 * c)) .* dot(((Kxx) .* zz - (Kzz) .* xx), S);
% mz = (1 ./ (2 * c)) .* dot(((Kyy) .* xx - (Kxx) .* yy), S);

coefficient_M = ((2 .* w.^4) / (3 .* c^3));
Imx = coefficient_M .* (abs(mx)).^2; % %磁偶极矩辐射能量
Imy = coefficient_M .* (abs(my)).^2;
Imz = coefficient_M .* (abs(mz)).^2;
Im = sqrt(abs(Imx).^2 + abs(Imy).^2 + abs(Imz).^2);
clearvars mx my mz cofficient_M Imx Imy Imz

%% 环形偶极矩
rr = sqrt(x .* x + y .* y + z .* z); % %位矢的模
rr = repmat(rr, 1, n); % %rr重复扩充n列,[m]
Tx = (1 ./ (10 * c)) .* dot((xx .* (xx .* Kxx + yy .* Kyy + zz .* Kzz) - 2 .* rr.^2 .* Kxx), S);
Ty = (1 ./ (10 * c)) .* dot((yy .* (xx .* Kxx + yy .* Kyy + zz .* Kzz) - 2 .* rr.^2 .* Kyy), S);
Tz = (1 ./ (10 * c)) .* dot((zz .* (xx .* Kxx + yy .* Kyy + zz .* Kzz) - 2 .* rr.^2 .* Kzz), S);
% Tx = (1 ./ (10 * c)) .* dot((xx .* (xx .* Kxx) - 2 .* rr.^2 .* Kxx), S);
% Ty = (1 ./ (10 * c)) .* dot((yy .* (yy .* Kyy) - 2 .* rr.^2 .* Kyy), S);
% Tz = (1 ./ (10 * c)) .* dot((zz .* (zz .* Kzz) - 2 .* rr.^2 .* Kzz), S);

coefficient_T = (2 .* w.^6) ./ (3 .* c^5);
ITx = coefficient_T .* (abs(Tx)).^2; % %环形偶极矩辐射能量
ITy = coefficient_T .* (abs(Ty)).^2;
ITz = coefficient_T .* (abs(Tz)).^2;
IT = sqrt(abs(ITx).^2 + abs(ITy).^2 + abs(ITz).^2);
clearvars rr Tx Ty Tz coefficient_T ITx ITy ITz

%% 电四极矩
Qxx = (1 ./ (i .* w)) .* dot(((xx .* Kxx) + (xx .* Kxx) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qxy = (1 ./ (i .* w)) .* dot(((xx .* Kyy) + (yy .* Kxx) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qxz = (1 ./ (i .* w)) .* dot(((xx .* Kzz) + (zz .* Kxx) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qyx = (1 ./ (i .* w)) .* dot(((yy .* Kxx) + (xx .* Kyy) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qyy = (1 ./ (i .* w)) .* dot(((yy .* Kyy) + (yy .* Kyy) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qyz = (1 ./ (i .* w)) .* dot(((yy .* Kzz) + (zz .* Kyy) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qzx = (1 ./ (i .* w)) .* dot(((zz .* Kxx) + (xx .* Kzz) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qzy = (1 ./ (i .* w)) .* dot(((zz .* Kyy) + (yy .* Kzz) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);
Qzz = (1 ./ (i .* w)) .* dot(((zz .* Kzz) + (zz .* Kzz) - 2/3 .* (xx .* Kxx + yy .* Kyy + zz .* Kzz)), S);

coefficient_Qe = ((w.^6) / (5 .* c^5));
IQxx = coefficient_Qe .* Qxx .* Qxx; % %电四极矩辐射能量
IQxy = coefficient_Qe .* Qxy .* Qxy;
IQxz = coefficient_Qe .* Qxz .* Qxz;
IQyx = coefficient_Qe .* Qyx .* Qyx;
IQyy = coefficient_Qe .* Qyy .* Qyy;
IQyz = coefficient_Qe .* Qyz .* Qyz;
IQzx = coefficient_Qe .* Qzx .* Qzx;
IQzy = coefficient_Qe .* Qzy .* Qzy;
IQzz = coefficient_Qe .* Qzz .* Qzz;
IQ = sqrt(abs(IQxx).^2 + abs(IQxy).^2 + abs(IQxz).^2 + abs(IQyx).^2 + abs(IQyy).^2 + abs(IQyz).^2 + abs(IQzx).^2 + abs(IQzy).^2 + abs(IQzz).^2);

%% 磁四极矩
RKx = conj(Kzz) .* conj(yy) - conj(Kyy) .* conj(zz); % %Rx(x,y,z)叉乘K(xx,yy,zz),x项;
RKy = conj(Kxx) .* conj(zz) - conj(Kzz) .* conj(xx); % %Rx(x,y,z)叉乘K(xx,yy,zz),y项;
RKz = conj(Kyy) .* conj(xx) - conj(Kxx) .* conj(yy); % %Rx(x,y,z)叉乘K(xx,yy,zz),z项;

Mxx = 1 ./ (3 .* c) .* dot((RKx .* xx) + (RKx .* xx), S);
Mxy = 1 ./ (3 .* c) .* dot((RKx .* yy) + (RKy .* xx), S);
Mxz = 1 ./ (3 .* c) .* dot((RKx .* zz) + (RKz .* xx), S);
Myx = 1 ./ (3 .* c) .* dot((RKy .* xx) + (RKx .* yy), S);
Myy = 1 ./ (3 .* c) .* dot((RKy .* yy) + (RKy .* yy), S);
Myz = 1 ./ (3 .* c) .* dot((RKy .* zz) + (RKz .* yy), S);
Mzx = 1 ./ (3 .* c) .* dot((RKz .* xx) + (RKx .* zz), S);
Mzy = 1 ./ (3 .* c) .* dot((RKz .* yy) + (RKy .* zz), S);
Mzz = 1 ./ (3 .* c) .* dot((RKz .* zz) + (RKz .* zz), S);

coefficient_Qm = ((w.^6) / (20 .* c^5));
IMxx = coefficient_Qm .* Mxx .* Mxx; % %磁四极矩辐射能量
IMxy = coefficient_Qm .* Mxy .* Mxy;
IMxz = coefficient_Qm .* Mxz .* Mxz;
IMyx = coefficient_Qm .* Myx .* Myx;
IMyy = coefficient_Qm .* Myy .* Myy;
IMyz = coefficient_Qm .* Myz .* Myz;
IMzx = coefficient_Qm .* Mzx .* Mzx;
IMzy = coefficient_Qm .* Mzy .* Mzy;
IMzz = coefficient_Qm .* Mzz .* Mzz;
IM = sqrt(abs(IMxx).^2 + abs(IMxy).^2 + abs(IMxz).^2 ++ abs(IMyx).^2 + abs(IMyy).^2 + abs(IMyz).^2 + abs(IMzx).^2 + abs(IMzy).^2 + abs(IMzz).^2);

%% 绘图 %%总量
figure('Renderer', 'painters', 'Position', [1041 304 420 * 1.65 420 * 1.65]);
plot(f / (1e12), Ip, '-o', 'Color', [0.3 0 0.5], 'LineWidth', 3.5, 'MarkerSize', 10, 'MarkerEdgeColor', [0.3 0 0.5], 'MarkerFaceColor', [0.49 1 0.63])
hold on
plot(f / (1e12), Im, '-d', 'Color', [0.5 0.5 0.5], 'LineWidth', 3.5, 'MarkerSize', 10, 'MarkerEdgeColor', [0.5 0.5 0.5], 'MarkerFaceColor', [0.49 1 0.63])
hold on
plot(f / (1e12), IT, '-v', 'Color', [0.5 0 0.3], 'LineWidth', 3.5, 'MarkerSize', 10, 'MarkerEdgeColor', [0.5 0 0.3], 'MarkerFaceColor', [0.49 1 0.63])
hold on
plot(f / (1e12), IQ, '-s', 'Color', [1 0 0], 'LineWidth', 3.5, 'MarkerSize', 10, 'MarkerEdgeColor', [1 0 0], 'MarkerFaceColor', [0.49 1 0.63])
hold on
plot(f / (1e12), IM, '-^b', 'Color', [1 0.5 0], 'LineWidth', 3.5, 'MarkerSize', 10, 'MarkerEdgeColor', [1 0.5 0], 'MarkerFaceColor', [0.49 1 0.63])
hold on
set(gca, 'FontSize', 28, 'Fontname', 'Helvetica', 'FontWeight', 'bold', 'linewidth', 3);
xlabel('Frequency (THz)', 'FontSize', 28, 'FontWeight', 'bold')
ylabel('Radiation   power(a.u.)', 'FontSize', 28, 'FontWeight', 'bold')
% set(gca,'YScale','log');% Y轴取对数显示
legend('\fontsize{28}\fontname{Helvetica}\itP', '\fontsize{28}\fontname{Helvetica}\itM', '\fontsize{28}\fontname{Helvetica}\itT', '\fontsize{28}\fontname{Helvetica}\itQ_e', '\fontsize{28}\fontname{Helvetica}\itQ_m')
legend('boxoff')
set(legend, 'Position', [0.721085860689599 0.520382406394848 0.155844153392883 0.388167376857127]);
box on

YW_Plot_Fast({f / (1e12), Ip, f / (1e12), Im, f / (1e12), IT, f / (1e12), IQ, f / (1e12), IM}, {'-o', '-d', '-v', '-s', '-^b'}, ["\fontsize{28}\fontname{Helvetica}\itP", '\fontsize{28}\fontname{Helvetica}\itM', '\fontsize{28}\fontname{Helvetica}\itT', '\fontsize{28}\fontname{Helvetica}\itQ_e', '\fontsize{28}\fontname{Helvetica}\itQ_m'], ["Frequency (THz)", 'Radiation   power(a.u.)'])
set(gca, 'yscale', 'log');
saveRadiativePower = [f' Ip' Im' IT' IQ' IM'];
